clc; clear; close all;

alpha = 1; sigma = 4; theta = 3.8;

I = 100;
 
%%%%%%%%%%%%%%%%%% Trade Fundamentals %%%%%%%%%%%%%%%%%%%%%%
b = 1*ones(I,1);
tau =  1*ones(I,I);
f =  0.0005*ones(I,I);
for i = 1:I
    tau(i,i) = 1;
    f(i,i) = 0.0005;
end
L = 10*ones(I,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

GG = 200;
grid = exp(linspace(log(40),log(1),GG));


for g = 1:GG

tau = grid(g)*ones(I,I);
for i = 1:I
   tau(i,i) = 1; 
end
    
w = ones(I,1);
y = ones(I,1);
P = ones(I,1);
 
tol = 1; iter = 0;

while tol > 10^(-5)||iter < 10
    
iter = iter + 1;

y = sigma*theta./(sigma*theta-(sigma-1)).*w;

onerho = repmat(b,1,I).^(theta).*(sigma.*repmat(w',I,1).*f).^(theta/(1-sigma)).*(sigma./(sigma-1).*repmat(w,1,I).*tau).^(-theta).*repmat(P',I,1).^theta.*repmat((L.*y)',I,1).^(theta./(sigma-1));
onerho = min(max(onerho,0),1);
star = (1-onerho);

lambda = repmat(L.*b.^theta,1,I).*(repmat(w,1,I).*tau).^(-theta).*f.^(1-theta./(sigma-1));
lambda = lambda./repmat(sum(lambda,1),I,1);

A = theta./(theta-(sigma-1)).*(sigma./(sigma-1)).^(-sigma).*(theta./(sigma*theta-(sigma-1))).^(theta./(sigma-1)-1).*L.^((theta-(sigma-1))./(sigma-1));
A = A.^(-1/theta);
    
P = A.*(sum(repmat(L.*b.^theta,1,I).*(repmat(w,1,I).*tau).^(-theta).*f.^(1-theta./(sigma-1)),1).^(-1/theta))';

wnew = sum(lambda.*repmat((L.*w)',I,1),2)./L;
wnew = wnew./wnew(1);

tol = sum(abs(wnew-w));
w = 0.8*w + 0.2*wnew;

end

G(:,g) = alpha.*(sigma-1)./(sigma*theta) - 2.*alpha.*(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(4*theta-2*(sigma-1))...
    .* sum(lambda.*onerho.*repmat((L.*y)',I,1)./repmat(L.*y,1,I),2);

s = 0.95;
Lorenz5(:,g) = (1-alpha.*(sigma-1)./(sigma*theta)).*s + alpha.*(theta-(sigma-1))./(sigma*theta).*sum((ones(I,I).*(s>star))...
    .*lambda.*repmat((L.*y)',I,1)./repmat(L.*y,1,I).*(theta./(theta-(sigma-1)).*(1-((1-s)./onerho).^((1-sigma)/theta+1))...
    -1 + (1-s)./onerho),2);
Share5(:,g) = 1 - Lorenz5(:,g);


s = 0.99;
Lorenz1(:,g) = (1-alpha.*(sigma-1)./(sigma*theta)).*s + alpha.*(theta-(sigma-1))./(sigma*theta).*sum((ones(I,I).*(s>star))...
    .*lambda.*repmat((L.*y)',I,1)./repmat(L.*y,1,I).*(theta./(theta-(sigma-1)).*(1-((1-s)./onerho).^((1-sigma)/theta+1))...
    -1 + (1-s)./onerho),2);
Share1(:,g) = 1 - Lorenz1(:,g);

rhostar(:,g) = onerho(1,2);
lambdastar(:,g) = 1-lambda(1,1);
end

plot(rhostar,G(1,:)./G(1,1)); hold on;
plot(rhostar,Share5(1,:)./Share5(1,1)); 
plot(rhostar,Share1(1,:)./Share1(1,1)); 
 
 
 
 